/*=========================================================================
 *
 *  Copyright NumFOCUS
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         https://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/
#ifndef itkHoughTransform2DLinesImageFilter_hxx
#define itkHoughTransform2DLinesImageFilter_hxx

#include "itkImageRegionIteratorWithIndex.h"
#include "itkDiscreteGaussianImageFilter.h"
#include "itkMinimumMaximumImageCalculator.h"
#include "itkCastImageFilter.h"
#include "itkMath.h"

namespace itk
{

template <typename TInputPixelType, typename TOutputPixelType>
HoughTransform2DLinesImageFilter<TInputPixelType, TOutputPixelType>::HoughTransform2DLinesImageFilter()
  : m_LinesList()

{}


template <typename TInputPixelType, typename TOutputPixelType>
void
HoughTransform2DLinesImageFilter<TInputPixelType, TOutputPixelType>::EnlargeOutputRequestedRegion(DataObject * output)
{
  // Call the superclass' implementation of this method
  Superclass::EnlargeOutputRequestedRegion(output);

  output->SetRequestedRegionToLargestPossibleRegion();
}


template <typename TInputPixelType, typename TOutputPixelType>
void
HoughTransform2DLinesImageFilter<TInputPixelType, TOutputPixelType>::GenerateOutputInformation()
{
  // Call the superclass' implementation of this method
  Superclass::GenerateOutputInformation();

  // Get pointers to the input and output
  const InputImageConstPointer input = this->GetInput();
  const OutputImagePointer     output = this->GetOutput();

  if (!input || !output)
  {
    return;
  }

  // Compute the size of the output image
  Size<2> size;

  size[0] = (SizeValueType)(std::sqrt(m_AngleResolution * m_AngleResolution +
                                      input->GetLargestPossibleRegion().GetSize()[0] *
                                        input->GetLargestPossibleRegion().GetSize()[0]));
  size[1] = (SizeValueType)m_AngleResolution;
  const typename InputImageType::RegionType region(input->GetLargestPossibleRegion().GetIndex(), size);

  output->SetLargestPossibleRegion(region);
}


template <typename TInputPixelType, typename TOutputPixelType>
void
HoughTransform2DLinesImageFilter<TInputPixelType, TOutputPixelType>::GenerateInputRequestedRegion()
{
  Superclass::GenerateInputRequestedRegion();
  if (this->GetInput())
  {
    const InputImagePointer image = const_cast<InputImageType *>(this->GetInput());
    image->SetRequestedRegionToLargestPossibleRegion();
  }
}


template <typename TInputPixelType, typename TOutputPixelType>
void
HoughTransform2DLinesImageFilter<TInputPixelType, TOutputPixelType>::GenerateData()
{
  // Get the input and output pointers
  const InputImageConstPointer inputImage = this->GetInput(0);
  const OutputImagePointer     outputImage = this->GetOutput(0);

  // Allocate the output
  this->AllocateOutputs();
  outputImage->FillBuffer(0);

  const double nPI = 4.0 * std::atan(1.0);

  ImageRegionConstIteratorWithIndex<InputImageType> image_it(inputImage, inputImage->GetRequestedRegion());

  Index<2> index;

  while (!image_it.IsAtEnd())
  {
    if (image_it.Get() > m_Threshold)
    {
      for (double angle = -nPI; angle < nPI; angle += nPI / m_AngleResolution)
      {
        index[0] =
          // m_R
          (IndexValueType)(image_it.GetIndex()[0] * std::cos(angle) + image_it.GetIndex()[1] * std::sin(angle));
        // m_Theta
        index[1] = (IndexValueType)((m_AngleResolution / 2) + m_AngleResolution * angle / (2 * nPI));

        if (index[0] > 0 && index[0] <= (IndexValueType)outputImage->GetBufferedRegion().GetSize()[0])
        // The preceding "if" should be replaceable with "if (
        // outputImage->GetBufferedRegion().IsInside(index) )" but
        // the algorithm fails if it is
        {
          outputImage->SetPixel(index, outputImage->GetPixel(index) + 1);
        }
      }
    }
    ++image_it;
  }
}


template <typename TInputPixelType, typename TOutputPixelType>
void
HoughTransform2DLinesImageFilter<TInputPixelType, TOutputPixelType>::Simplify()
{
  // Get the input and output pointers.
  const InputImageConstPointer inputImage = this->GetInput(0);
  const OutputImagePointer     outputImage = this->GetOutput(0);

  if (!inputImage || !outputImage)
  {
    itkExceptionMacro("Update() must be called before Simplify().");
  }

  // Allocate the simplify accumulator.
  m_SimplifyAccumulator = OutputImageType::New();
  m_SimplifyAccumulator->SetRegions(outputImage->GetLargestPossibleRegion());
  m_SimplifyAccumulator->SetOrigin(inputImage->GetOrigin());
  m_SimplifyAccumulator->SetSpacing(inputImage->GetSpacing());
  m_SimplifyAccumulator->SetDirection(inputImage->GetDirection());
  m_SimplifyAccumulator->AllocateInitialized();

  Index<2> index;
  Index<2> maxIndex;

  typename OutputImageType::PixelType value;
  typename OutputImageType::PixelType valuemax;

  ImageRegionConstIteratorWithIndex<InputImageType> image_it(inputImage, inputImage->GetRequestedRegion());

  const double nPI = 4.0 * std::atan(1.0);

  while (!image_it.IsAtEnd())
  {
    if (image_it.Get() > m_Threshold)
    {
      // Look for maximum along the curve and remove the curve at the same time.
      valuemax = -1;
      maxIndex[0] = 0;
      maxIndex[1] = 0;
      // Compute the distance from the corner to the line (first dimension),
      // and the angle between the X axis and the normal to the line (second
      // dimension).
      for (double angle = -nPI; angle < nPI; angle += nPI / m_AngleResolution)
      {
        // m_R
        index[0] =
          (IndexValueType)(image_it.GetIndex()[0] * std::cos(angle) + image_it.GetIndex()[1] * std::sin(angle));
        // m_Theta
        index[1] = (IndexValueType)((m_AngleResolution / 2) + m_AngleResolution * angle / (2 * nPI));

        if (outputImage->GetBufferedRegion().IsInside(index))
        {
          value = outputImage->GetPixel(index);
          if (value > valuemax)
          {
            valuemax = value;
            maxIndex = index;
          }
        }
      }
      m_SimplifyAccumulator->SetPixel(maxIndex, m_SimplifyAccumulator->GetPixel(maxIndex) + 1);
    }
    ++image_it;
  }

  ImageRegionConstIteratorWithIndex<OutputImageType> accusimple_it(m_SimplifyAccumulator,
                                                                   m_SimplifyAccumulator->GetRequestedRegion());
  ImageRegionIteratorWithIndex<OutputImageType>      accu_it(outputImage, outputImage->GetRequestedRegion());

  while (!accusimple_it.IsAtEnd())
  {
    accu_it.Set(accusimple_it.Get());
    ++accu_it;
    ++accusimple_it;
  }
}


template <typename TInputPixelType, typename TOutputPixelType>
auto
HoughTransform2DLinesImageFilter<TInputPixelType, TOutputPixelType>::GetLines() -> LinesListType &
{
  // If the filter has not been updated.
  if (this->GetMTime() == m_OldModifiedTime)
  {
    return m_LinesList;
  }

  m_LinesList.clear();

  if (m_NumberOfLines > 0)
  {
    // Blur the accumulator in order to find the maximum.
    using InternalImagePixelType = float;
    using InternalImageType = Image<InternalImagePixelType, 2>;

    const OutputImagePointer outputImage = this->GetOutput(0);

    if (!outputImage)
    {
      itkExceptionMacro("Update() must be called before GetLines().");
    }

    // Convert the accumulator output image type to internal image type.
    using CastImageFilterType = CastImageFilter<OutputImageType, InternalImageType>;

    const auto castImageFilter = CastImageFilterType::New();
    castImageFilter->SetInput(outputImage);

    using GaussianFilterType = DiscreteGaussianImageFilter<InternalImageType, InternalImageType>;
    const auto gaussianFilter = GaussianFilterType::New();

    // The output is the accumulator image.
    gaussianFilter->SetInput(castImageFilter->GetOutput());
    gaussianFilter->SetVariance(m_Variance);
    gaussianFilter->Update();
    const InternalImageType::Pointer postProcessImage = gaussianFilter->GetOutput();

    using MinMaxCalculatorType = MinimumMaximumImageCalculator<InternalImageType>;
    auto                                        minMaxCalculator = MinMaxCalculatorType::New();
    itk::ImageRegionIterator<InternalImageType> it_input(postProcessImage,
                                                         postProcessImage->GetLargestPossibleRegion());


    itk::Index<2> index;

    unsigned int lines = 0;

    // Find maxima
    do
    {
      minMaxCalculator->SetImage(postProcessImage);
      minMaxCalculator->ComputeMaximum();
      InternalImageType::PixelType max = minMaxCalculator->GetMaximum();

      if (max <= 0)
      {
        // When all pixel values in 'postProcessImage' are zero or less, no more lines
        // should be found. Note that a zero in 'postProcessImage' might correspond to a
        // zero in the accumulator image, or it might be that the pixel is within a
        // removed disc.
        break;
      }

      for (it_input.GoToBegin(); !it_input.IsAtEnd(); ++it_input)
      {
        if (Math::ExactlyEquals(it_input.Get(), max))
        {
          // Create the line.
          LineType::LinePointListType list; // Insert two points per line.

          const double radius = it_input.GetIndex()[0];
          const double teta = ((it_input.GetIndex()[1]) * 2 * Math::pi / this->GetAngleResolution()) - Math::pi;
          const double Vx = radius * std::cos(teta);
          const double Vy = radius * std::sin(teta);
          const double norm = std::sqrt(Vx * Vx + Vy * Vy);
          double       VxNorm = Vx / norm;
          double       VyNorm = Vy / norm;

          if (teta <= 0 || teta >= Math::pi / 2)
          {
            if (teta >= Math::pi / 2)
            {
              VyNorm = -VyNorm;
              VxNorm = -VxNorm;
            }

            LinePointType p;
            p.SetPositionInObjectSpace(Vx, Vy);
            list.push_back(p);
            p.SetPositionInObjectSpace(Vx - VyNorm * 5, Vy + VxNorm * 5);
            list.push_back(p);
          }
          else
          {
            LinePointType p;
            p.SetPositionInObjectSpace(Vx, Vy);
            list.push_back(p);
            p.SetPositionInObjectSpace(Vx - VyNorm * 5, Vy + VxNorm * 5);
            list.push_back(p);
          }

          // Create a Line Spatial Object.
          const LinePointer line = LineType::New();
          line->SetId(lines);
          line->SetPoints(list);
          line->Update();

          m_LinesList.push_back(line);

          // Remove a black disc from the hough space domain.
          for (double angle = 0; angle <= 2 * Math::pi; angle += Math::pi / 1000)
          {
            for (double length = 0; length < m_DiscRadius; length += 1)
            {
              index[0] = (IndexValueType)(it_input.GetIndex()[0] + length * std::cos(angle));
              index[1] = (IndexValueType)(it_input.GetIndex()[1] + length * std::sin(angle));
              if (postProcessImage->GetBufferedRegion().IsInside(index))
              {
                postProcessImage->SetPixel(index, 0);
              }
            }
          }
          minMaxCalculator->SetImage(postProcessImage);
          minMaxCalculator->ComputeMaximum();
          max = minMaxCalculator->GetMaximum();

          ++lines;
          if (lines == m_NumberOfLines)
          {
            break;
          }
        }
      }
    } while (lines < m_NumberOfLines);
  }

  m_OldModifiedTime = this->GetMTime();
  return m_LinesList;
}


template <typename TInputPixelType, typename TOutputPixelType>
void
HoughTransform2DLinesImageFilter<TInputPixelType, TOutputPixelType>::PrintSelf(std::ostream & os, Indent indent) const
{
  Superclass::PrintSelf(os, indent);

  os << indent << "Threshold: " << m_Threshold << std::endl;
  os << indent << "Angle Resolution: " << m_AngleResolution << std::endl;
  os << indent << "Number Of Lines: " << m_NumberOfLines << std::endl;
  os << indent << "Disc Radius: " << m_DiscRadius << std::endl;
  os << indent << "Accumulator blur variance: " << m_Variance << std::endl;
  itkPrintSelfObjectMacro(SimplifyAccumulator);

  os << indent << "LinesList: " << std::endl;
  unsigned int i = 0;
  auto         it = m_LinesList.begin();
  while (it != m_LinesList.end())
  {
    os << indent << '[' << i << "]: " << *it << std::endl;
    ++it;
    ++i;
  }

  os << indent << "OldModifiedTime: " << NumericTraits<ModifiedTimeType>::PrintType(m_OldModifiedTime) << std::endl;
}
} // end namespace itk

#endif
